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ABSTRACT 


An algorithm for the computer-aided design of passive electrical 
filters using minimization techniques is presented. The configuration, 
or topology of the filter and its desired frequency response are speci- 
fied. A measure of the deviation between the desired and the actual 
frequency response is minimized using the gradient projection search 
method. Examples are presented which demonstrate the performance of 


the minimization procedure. 
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IT. INTRODUCTION 


The development of extremely high-speed computers and efficient 
algorithms has opened the door to the use of iteration techniques that 
heretofore have been prohibitively costly in time and/or money. This 
high-speed capability makes it possible to consider a computer-aided 
solution to the problem of designing an electrical filter to match — 
as closely as possible —a given frequency response, wherein many 
hundreds ~— sometimes thousands — of successive solutions of the network 
frequency response are required. Use of the designing engineer's expe- 
rience in adaptation of filters whose response is already known, but is 
not quite that desired, usually demands laborious hand calculation. In 
addition, there is occasional difficulty in achieving the required 
“closeness of fit" to the desired response. Some frequency response 
requirements may lie beyond the capability of synthesis techniques, 
being solvable only by experience coupled with trial and error proce- 
dures — which may never produce the desired accuracy. 

With the availability of a computer program capable of systematically 
minimizing deviations from the required response, design time can be cut 
dramatically. The ease with which this tool can produce results has 
multiple benefits: 


1. The engineer is free to concentrate on design without the tedium 
of repetitious computation. 


2. The designer can investigate numerous circuit configurations to 
determine one which "best" fits the needs of simplicity, cost, and 
quality of performance. 


3. The computer's speed enables the designer to quickly obtain an 
answer. 


11. OPTIMIZATION THEORY 


A. MINIMIZING THROUGH OPTIMIZATION 

Optimization techniques are used to provide the means of systemat- 
ically and effectively developing the soiution to the probiem. Classical 
methods for designing filters to meet a given frequency response are 
highly time-consuming, and usually produce only approximations of the 
desired response, perhaps matching the general shape, and a specified 
bandwidth [Ref. 1]. Optimization techniques utilize the designer's 
particular definition for "closeness of fit" of the designed to the 
specified frequency response and then proceed through successive adjust- 
ments of the variable parameters to minimize the deviation. The function 
defined by the designer to specify this closeness of fit will be referred 
to as "the cost function". 

A specific filter configuration is required by tne optimization 
process; no provision is made to permit computer-determined alteration of 
the circuit in an effort to improve the overall closeness of fit. This 
decision is left to the engineer to be carried out during his interaction 


with the program. 


B. MINIMIZATION TECHNIQUES 

Three methods for minimization were considered: direct, or pattern 
search; the gradient method; and second-order methods. 

1. Pattern Search 

The pattern-search method makes use of successive perturbations 

of the variable parameters to determine when a lower value of the cost 
function is achieved. When a reduction in value of the cost function is 
determined, the new point is then used to again perturb the parameters, 
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taking larger steps when continued successes are achieved. If a failure 
to decrease the cost function occurs, the step size is reduced until 
success is achieved. This is the simplest of the three methods, re- 
quiring only values of the cost function, but is somewhat sensitive to 
starting point when local minima are present [Ref. 2]. 
2. Gradient Method 

The gradient method makes use of the negative of the first 
partial derivative of the cost function with respect to each of the 
variable parameters; this determines the direction of steepest descent. 
A step of specified distance is taken in the steepest-descent direction 
to decrease the vaiue of the cost function. 

A closely related method, that of gradient projection [Ref. 3], 
makes use of the negative gradient at the starting point to extend a 
vector in this "steepest-descent" direction until one of the parameter 
constraints is reached. The smallest cost function value along this 
extended vector is then determined. The negative gradient vector is 
computed for this new point, and the previous step is repeated until the 
minimum is reached. The algorithm will always move in the negative 
gradient direction unless coing so causes a violation of the constraints, 
in which case the negative gradient vector is projected onto the con- 
straint boundary, and the next step is taken in the direction of this 
projection. 

3. Second-Order Methods 

Second-order methods involve the use of both the first and 
second partial derivatives of the cost function, and were not considered 
once the difficulty in generating the second partial derivative for 


filter networks was recognized. 


C. PROBLEM CHARACTERISTICS 
1. The Cost Function 
The general format for the cost function, J(p), used in the 


optimization scheme is defined as: 


Jp) = s¢ hlw) fe(p)|* du - (2.1) 


Onin 


Here p is the vector of variable circuit parameters; h(w) is a weighting 
function used to force the solution to more closely fit the prescribed 
response in certain critical areas; e(p) is the error between the desired 
and the actual frequency response; N represents the power to which e(p) 
is raised to determine the type of closeness of fit (e.g., N = 2 for mean- 


square), and normally takes on the value of 1, 2 or 43 w. 


naan and a are 


ax 
the lower and upper frequencies over which the integration is performed. 


The integration must be written in discrete form for use in a 


computer program; thus the following approximation is used; 


+kAw)-R. (wt +kAw) | ee 


Dwg tkaw) [Ral “min 


Omin a 


In this equation Aw is the difference between successive frequency points; 
K represents the number of frequency points to be matched; N is as 


described above, and h(w_- +kaw) is the discrete equivalent of h(w); 


min 
R lw. tkaw) is the desired frequency response at the ge frequency 


point, while Ro (us tkaw) is the actual response. It is this function, 


J(p), that the optimization process is to minimize. If a perfect fit to 


~ 


the desired frequency response is achieved, the value of J(p) will be 


~ 


zero, its absolute minimum. 


a Merininatrony unl cer ta 
Iteration of an optimization program can be terminated i% several 
ways. Three of these are: 
1. When further iteration produces no meaningful improvement in close- 
ness of fit to the desired frequency response, iteration is terminated. 


A suitable equation is: 


Oe Ne (2.3) 


_ ~ 


where e« is some relatively small positive number. If ce is selected to be 
too small, excessive iterations will usually be required to produce two 
successive values of J(p) close enough together to terminate the procedure, 
resulting in wasted computer time. If a large value is chosen for e« , the 
program may terminate far from a minimum. 

2. When the norm of the gradient vector is quite close to zero, either 
a minimum, a maximum or a saddle point has been found; in a constrained 
problem, maximums are less likely. Numerous problems were run, and 


neither a maximum nor a saddle point ever resuited as an optimization 


solution. This criterion can be written: 


where y is a suitably small positive number. 
3. A third criterion, not written into the computer program used in 
this investigation, but one which might prove useful to the designer in 


utilization of the program is: 


, 


where 6 1S some suitably small positive number. Inclusion of this 
inequality into the computer program would cause termination once the 
cost function reaches a value less than 6. 
3, Constraints 
In many practical design problems, there is some idea of a 
reasonable range of constraints for each variable parameter, therefore 
it is assumed that constraints of the form shown below can be imposed on 


the n variable parameters: 


where a = (a) sAgo-ee sa) represents the vector lower limit, and 


~s 


b = (b) sb5.-++ sb) is the vector upper limit. In the case of several 


variables, the constraints influence the rate of convergence, especially 
for the gradient-projection method. For example, a reduction of con- 

Suro cs, + POltaUr <P, <'. to 0. <p, < -5 , where n = 5, will reduce 
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a five-dimensional space by 96.9 percent (using (1/2)~ = .03124), leaving 


a smaller region in which to search for the minimum. 


le 


III. THE OPTIMIZATION PROGRAM 


A. THE FREQUENCY ANALYSIS PROGRAM 
Several existing computer programs for calculating the frequency 
response of a given circuit were studied. Some of these were: 


Calahan 
Sceptre 
Ecap 
Cornap 


The Calahan Network Analysis Program was chosen as being best able to 
provide the required data, in the desired form, with high accuracy and 
reasonable speed. 
1. The Basic Program 

The Linear Network Analysis Program by D. A. Calahan [Ref. 4] is 
a large general purpose program capable of producing a variety of circuit 
information. Inputs to the program are: either the network transfer 
function or a topological description; keys to indicate the desired output; 
high- and low-frequency limits and the number of points to be determined, 
if the frequency response is desired; initial and final values, initial 
and final time, and the time increment if a time response is desired. 
Provisions are made for varying one element for as many as eight differ- 
ent values. The output from the program includes the coefficients of the 
network function, the poles and zeros, the frequency response, the time 
response (to an arbitrary input), and the symbolic form of the network 
function. The program is capable of accepting circuits with as many as 
thirty nodes, one hundred passive elements, and twenty active elements. 


Printer-piot graphs of frequency and time response are produced when 


TS 


desired. Negative element vaiues, as wel! as negative frequencies, are 
acceptable to the program. 
2. Accuracy 

Examination of the accuracy of the frequency response produced 
by the Caiahan analysis was accomplished in the following fashion. 
Specific circuits whose frequency response could be readily calculated 
using formulae (Butterworth and Chebychev low-pass filters in particular 
were used, producing the expected frequency response using formulae 
found in [Ref. 5]) were tested using the nominai element values. The 
Calahan results compared favorably with the answers from the formulae. 
For example, analysis of a Tive-element Butterworth filter gave the 


following results: 


Formula Calanan 

Freq(Hz ) Gain(db) Gain(db) 
el -2.10196 -2.1018) 
16 =Sehuoe/ -3,10314 

es PT) -4 30471] ~4. 3046} 
ahs -5.65472 -5 65465 
19 -/ 09729 =7.09723 
20 -8 5844) -8 . 58439 


3. Modifications 
Because only the frequency response was desired as an output from 
Calahan, and this on a highly repetitive basis, modifications to the 
program were made. Several subroutines and the main program were removed; 
a new main program was substituted so that the solution to the minimization 
problem could be generated by repeated use of the Calahan subroutines to 
determine the required frequency response. Two of the remaining Calahan 


subroutines were modified to conform to the desired iteration scheme. 


B. THE MINIMIZATION PROGRAM 


The gradient, or steepest~descent method, was chosen to be used for 
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minimization of the previously specified cost function. (The pattern 
search solution to this problem is the subject of a Naval Postgraduate 
School thesis by Captain James Lau, USMC. Comparative results are 
discussed in Chapter VI.) 
1. Determination of the Gradient Vector 
Some thought was given to the best method of producing the vector 


3J(p) 


xD - The literal matrix of partial derivatives was manually computed 


for three- and five-element filters, and it was quickly seen that this 





method is much too cumbersome and time-consuming, and will greatly reduce 
the ease with which the designer can use the program. As an example, 
given below is one partial for a three-element filter (two variable 
elements, L and C3; R=1). 

K 
» 


20 LC“(kaw) “+ {kaw) “-c( kaw)” J 


at ee (kaw) + (°-21.0) (kaw) “09 


oJ (LC) 
OL 





1 


—% a. . a Lo - —-——7} - (3.1) 


°C? (kaw) #(Lo-2LC)(kaw)“+1 = 1#( Kaw) 


It was then decided to investigate the method of finite-difference 
approximation. This method makes use of the definition of the partial 


derivative, 





ad (p) 7 J(Py+APy sPosece oP.) ~ JP] sPyo--+2P,) 3.2) 
ee Se ee : 
sail Ap,>0 =] 


The finite-difference approximation assigns a smali but finite value to 


AP, and determines the approximate partial derivative as: 


aJ(p) (Py AD, sPys--+9P,) - I(Py oPy 90-0 sP, ) 


~w 


ap AP; 
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Experimentation with various networks, using the Calahan subrou- 
tines to produce the frequency response indicated that a Ap of 0.0] 
percent gives reasonably accurate results. This accuracy was determined 
in two ways: by comparison with the values of the partial derivative 
vector for a two-variable-element filter produced by literal equations, 
one of which is equation (3.1); by use of discrete approximation to 
produce the gradient vector at the known absolute minimum of various 
larger networks. The first method produced normalized discrete approx- 
imation vectors within 1.5 percent of the true partial derivative vector 
in each of the ten points compared. The second method utilized larger 
networks whose actual frequency response was known, as well as the element 
values necessary to produce this response. When the correct element 
values were entered into the optimization program, each element of the 
gradient vector was less than 5 x ie. in all examples tested. 

2. The Gradient Method 

Initially, a number of circuit optimization problems were tried, 

using the gradient method. In this program, new values for the variables 


are produced using the following iteration equation: 


Each iteration produces a new set of parameter values, where the normal- 
ized gradient vector is "stepped" a distance g. When the new value for 
the cost function is greater than the previous value, g is halved and the 


step is repeated until a smaller function value is obtained. This method 
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produces good results on filters with up to five variable parameters, but 
required many iterations that reduce the value of the cost function very 
Slowly as the minimum is approached. A better method was therefore 
sought. 

2. The Gradient-Projection Method 

In hope of eliminating some of the disadvantages encountered in 
the gradient method, the gradient-projection method introduced by J. B. 
Rosen was investigated. This algorithm showed promise of more rapid 
convergence. Additionally, a program using this method was locally 
available. This is the Gradient Projection Program written by Dr. D. E. 
Kirk [Ref. 6]. The process is an iterative numerical method for locating 
the extremum of a function of many variables which are required to satisfy 
linear constraints. Unconstrained variables cannot be used with this 
method, as will be shown. 

A simple two-dimensional example will serve to illustrate the use 
of this method. Let f(p) be a differentiable function of two variables 
p, and p,; F (Pp, »P5) represents the value of the function at the point 
(P]»Po)- The objective is to find the point (P)*5po*) where f(p*) is a 


minimum, and which satisfies the constraints: 


0. <P; Saale 


Figure 1(a) represents a possible contour plot of f(p) as a function of 


~ 


P, and Po; the region which satisfies the constraints lies within the 


rectangle bounded by the lines Cl, C2, C3 and C4. The starting point is 


at p(O) = (2,1). At this point, the negative gradient vector, - aT) 


1s calculated, and extended until a constraint boundary is reached, in 


re, 





Figure | 


Contour Maps Demonstrating the Gradient Projection Method 


this case C4. The gradient at p(1), the intersection of C4 and the ex- 


tended vector 2) is calculated. The inner produced of - “Hes 
and - aore 1S negative, indicating that the minimum value of f(p) 


along this extended vector is somewhere between p(O0) and p(1)3; therefore 


interpolation is necessary. Successive interpolation locates the point 


p(2), where - af (2) 1s orthogonal to the originally extended vector. 


i 3p 
Orthogonality is indicated when the inner product of - 2f(C) iil _ af(2) 


dp dp 
1s zero. The negative gradient vector - ate) is then extended to the 


nearest constraint boundary, and the process continues until the minimum 
is reached at p(6), as shown in Figure l(a). 
Figure 1(b) shows a different hypothetical contour plot for f(p), 


we 


using the constraints: 


© 
ras 
Tw 
——_s 
fn 
a | 
(SE | 


Oo oe 


In this case, the absolute minimum lies outside the allowable constrained 
region; the minimum therefore lies along a boundary. This situation 
causes the constraint boundary C3 to be "active"; that is, it separates 
the acceptable region (that region within the constraint boundaries) from 


the absolute minimum of the function f(p). Once again, the starting point 


is p(O) = (2,1). As before, the inner product of - TAO) and - Jas 


is negative, indicating that interpolation is necessary. The point p(2) 


is found and its negative gradient extended to the boundary C3 where, 


unlike the preceeding step, the inner product of - ate) ane - wa! 1s 
positive, indicating that p(3) is the minimum along the extension of 


40) and that no interpolation is necessary. However, a change of p 


in the direction of - 2013) wourd cause a violation of the constraints; 


therefore p is changed aiong the vector represented by the projection of 


Re 





P A along the boundary C3. This projection 1s denoted by P[- 2113)), 
and is extended until another constraint is intersected, this time C4. 
Here the inner product of P[- ath) and - tA) negative, requiring 
interpolation. Successive interpolations locate p(5), where P[- heh 
and - 3f(5) are orthogonal. This is the minimum within the acceptable 


ap 
region. It should now be evident that unconstrained parameters cannot 


be used with this method. 

In Figures 1(a) and (b), the functions to be minimized are uni- 
modal. If, however, a given function is not unimodal within the acceptable 
region, then more than one minimum can exist, and the minimum found by the 
gradient-projection method may depend on the starting point. Figure 2 
demonstrates this, showing solutions for two different starting points 
(gradient vectors are not labeled in order to enhance clarity). A point 
of interest is demonstrated by the projection vector from the second 
starting point: the actual minimum along this projec*ion vector is 
located at A; however, the inner-product method indicates that no inter- 
polation is necessary, and accepts point p'(1) as the next starting point. 
This clearly demonstrates the requirement of unimodality for a function 
to be minimized by this interpolation method if the absolute minimum is 
desired, and multiple starting points are to be avoided. 

A more detailed explanation of this method, including the matrix 
equations required to Ws solved in the computer program, is contained in 
Reference 6. 

4. Cost-Function Complexity 

The inner-product interpolation scheme incorporated in the aradient 
projection program produces accurate results when a convex function is to 
be minimized; a convex function 1s unimodal. In using this interpolation 
method to optimize the frequency response of a filter network, it was found 
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that the scheme was not functioning so that ever-decreasing values of the 
cost function were obtained. This helped confirm what was already sus- 
pected: the cost function, J(p), as applied to filter network frequency 
response may not be a convex function. An actual contour map of J(p) as 
a function of two variable parameters is shown in Figures 3(a) and (b). 
This map was produced by solving for the cost function over a grid of 
values of the two variables and joining points of equal cost-function 
value. The map shows that while only one minimum exists, the function 
J(p) is obviously not convex. Higher-dimensional maps cannot be easily 
drawn, but experience in minimizing filters with several variable elements 
indicates that more than one minimum may be present. Therefore it is 
concluded that the cost function, J(p), as used in this procedure is not 
convex, and is possibly, in general, not unimodal. No conclusive results 
in this area were sought, as this was not the object of the investigation. 
5. The Function-Search Interpolation Procedure 

Armed with the knowledge of the cost function's complexity, a new 
interpolation method was sought. This method would ideally achieve the 
following goals: 

Goal 1 Continually decrease the value of the cost function, 
never accepting any new point unless it yields a lower value of J(p). 

Goal 2. Make efficient use of computer time. 

Goal 3. Increase the likelihood of converging to the absolute 
minimum, if local minima should be present 

Investigation of various possible methods produced one that showed 
excellent promise of fulfilling all three requirements. This interpolation 


does not use successive values of the gradient vector, but successive 


values of the cost function. A simple example is the easiest means of 
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that the scheme was not functioning so that ever-decreasing values of the 
cost function were obtained. This helped confirm what was already sus- 
pected: the cost function, J(p), as applied to filter network frequency 
response may not be a convex function. An actual contour map of J(p) as 
a function of two variable parameters is shown in Figures 3(a) and (b). 
This map was produced by solving for the cost function over a grid of 
values of the two variables and joining points of equal cost-function 
value. The map shows that while only one minimum exists, the function 
J(p) is obviously not convex. Higher-dimensional maps cannot be easily 
drawn, but experience in minimizing filters with several variable elements 
indicates that more than one minimum may be present. Therefore it is 
concluded that the cost function, J(p); as used in this procedure is not 
convex, and is possibly, in general, not unimodal. No conclusive results 
in this area were sought, as this was not the object of the investigation. 
5. The Function-Search Interpolation Procedure 

Armed with the knowledge of the cost function's complexity, a new 
interpolation method was sought. This method would ideally achieve the 
following goals: 

Goal 1. Continually decrease the value of the cost function, 
never accepting any new point unless it yields a lower value of J(p). 

Goal 2. Make efficient use of computer time. 

Goal 3. Increase the likelihood of converging to the absolute 
minimum, if local minima should be present 

Investigation of various possible methods produced one that showed 
excellent promise of fulfilling all three requirements. This rRenEoTEL Ten 


does not use successive values of the gradient vector, but successive 


values of the cost function. A simple example is the easiest means of 
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explaining the interpolation procedure. Using a filter with two variable 
elements (so that geometrical representations are possible), a contour 


map of J(p) as a function of the variables might be as shown in Figure 4. 
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Hypothetical Contour Map Showing the Segmented 
Extension of the Negative Gradient Vector 

At the base point where the cost is J(1), the negative gradient is calcu- 
lated and extended until it meets the first constraint boundary, p, (min), 
as would have been done with the inner product interpolation method. The 
length of the extended vector will be referred to as Tae Now, however, 
the line in the direction of the extended vector is subdivided into an 
even number of intervals (four intervals and five points will be used in 
this example). The cost function is determined for each of the remaining 
points, J(2) through J(5). The smallest value of J and its corresponding 
point number are determined. If J(5) is the smallest cost-function value, 
no interpolation is performed; if J(1) is the smallest, the interval 
between J(1) and J(2) is investigated, since the gradient vector's direc- 
tion indicates that some point with a smaller value than J(1) lies in this 


direction. If an interior point is smallest, interpolation must be © 
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accomplished. As is seen in Figure 4, J(2) has the smallest cost-function 
value. The interpolation scheme then passes a parabola through J(2) and 
its adjacent points (J(1) and J(3)), determines the minimum point of this 
parabola along the line in the direction of the extended vector, T, and 
evaluates the cost function, J(6), at this new point, as shown in Figure 
5(a). The purpose of approximating the actual curve of the cost function 
by a parabola, is that for reasonably well-behaved functions, a parabola 
joining three relatively close points on the curve will closely approxi- 
mate the actual curve. This helps to increase the speed with which the 
actual minimum is found. 

A detailed explanation of the procedure for finding the parabolic 
minimum is contained in Appendix A. 

Having solved for J(6) as shown in Figure 5(a), the points J(2), 
J(6), and J(3) are then used to solve for another parabolic minimum 
(Figure 5(b)). This process continues until the difference in extended 
vector length of the location of two successive minima is within five 


percent of ae In this case, location of parabolic minima would termi- 


x" 
nate if T, - Tp = 0.05 Lieven, The interpolation is then considered 
complete, and the lowest cost function value achieved is used as the new 
starting point. 

In numerous examples, this interpolation process has succeeded in 
quite reasonably accomplishing each of the desired goals. 

Goal 1. The requirement to continually decrease the value of the 
cost function was achieved. 

Goal 2. Computation time has been reduced over the inner- product 
interpolation method. Approximately the same number of points along the 


gradient vector are necessary in the function-search method, but only the 


value of the cost function, and not the gradient vector, is required at 


Zo 


Actual curve _— 
Parabolic estimate —— 
Actual minimum x 
Parabolic minimum o 


Distance Along the 
Extended Vector, T 





J(p) 
20 
a 9) = 
10 _ (6) —— J(3) 
J(7) 


T, Te T T, Distance Along the 
/ Extended Vector, T 


Figure 5 


Parabolic Estimation Diagrams 
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each point. This eliminates n(m - 2) cost-function calculations, where 

n is the number of variables, and m is the number of points considered 

along the projection vector (the function search method still requires 

the gradient vector at the first and last points, accounting for the -2). 
Goal 3. The likelihood of converging to the absolute minimum 

appears to have been increased by use of this method. Results from its 

use in circuit frequency-response minimization indicate that such may be 

the case. Appendix B cites two examples (many more have been noted) 

where a distinct "break" occurred in what appeared to be the final approach 

stages where the gradient procedure converges to a minimum point. A sudden 

and reasonably large change in parameter values and decrease in cost- 

function value take place, something that does not occur using gradient 


methods that do not make use of the extended vector. 


2/ 


WW. RESULTS 


A general result is first presented, followed by specific examples 
illustrating results for several of the design problems investigated. 
The final section of the chapter discusses a problem encountered during 


the investigation. 


A. GENERAL RESULTS 
Since only a limited number of filters were investigated, the results 
are not ‘conclusive but serve to demonstrate that optimization of the 


frequency response of a given filter is feasible using the techniques 


previously described. Although all the examples are of the passive variety, 


use of active elements should be feasible using models for the active 


elements required for the Calahan subroutines described in Reference 3. 


By SPRET RIC RESULTS 

In this section, specific results for those filters investigated that 
representatively demonstrate the capability of the optimization program 
are presented with circuit diagrams, program solutions and graphical 
results. 

1. Example 1 

Figure 6 contains the results of the optimization of a six-element 

low-pass filter. The desired curve is the Butterworth filter response 


feet . eae 
|Gain| = 20 Logi [(1 + oop les (4.1) 


where n is the number of reactive elements. In this case, n= 5. The 


actual response deviates from the desired only slightly. A constant 


28 


Circuit Configuration 





Starting Values of the Variable Parameters 
L. =e a | eei= 10 C, = 2.4 C, = 0.9 
Final Values of the Variable Parameters 


L, = qqglte5 57 Ly =e O9S L. = .500 C, = 1.478 C,, = 1.064 


Number of Iterations — 43 Computer Run-time — 136 seconds 


3) G40) -9 


Termination Caused by < 10 


Graphical Comparison of Desired and Actual Frequency Response 





Desired — 
Actual 


f (Hz ) 


Figure 6 


Data for Example | 
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weighting matrix (h. = 1., 1 = 1,...3K) was used in this example, and in 
all other examples unless otherwise stated. 

Figure 7 demonstrates the rapidity with which the program improves 
the actual response. The filter response is shown in comparison to the 
desired response after the indicated number of iterations for the problem 
described in Figure 6. After ten iterations the actual response is a 
relatively good approximation to the desired response. 

2. Example 2 

Figure 8 contains the results of the optimization of a five- 
element low-pass filter, using straight-line approximation for the desired 
frequency response. Since this is not a realizeable response, it is only 
possible to approximate it with the actual response. The results are 
reasonably close to the desired shape and bandwidth. If it were desired 
to more closely match the 3-db bandwidth, a weighting matrix could be 
used that places greater emphasis on the points where the desired response 
is -3 decibels. 

3, Eee 

Figure 9 contains the results of the optimization of an eight- 

element band-pass filter, wherein six elements were varied. The desired 


response was realizeable with the specified network configuration, and 


was very closely approximated, except in the frequency interval ((.14,0.22) 


where the deviation was noticeable. The bandwidth and general shape are 
reasonably accurate. This same filter circuit was optimized using a 
weighting matrix where both the upper and lower 3-db points were weighted 
with a value five times that of the weight of the other points. The devi- 


ation at the 3-db points was improved by a factor of about three. 
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Circuit Configuration 





Starting Values of the Variable Parameters 


= 2.0 Cae 60 Crue c.0 


Final Values of the Variable Parameters 


L, = 1.434 L,, = 1.093 C, = 1.635 Cy = 0.385 


Number of Iterations — 40 Computer Run-time — 132 seconds 


Termination Caused by gf) - g (41) < i0° 


Graphical Comparison of Desired and Actual Frequency Response 


|Gain| 
chy | # 
<a Desired — 
a Actual ---- 
5 Ns 
-10 
-15 
-20 £ (Hz) 
0.10 OMS 0226 O25 
Figure 8 


Data for Example 2 
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Circuit Configuration 
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Final Values of the Variable Parameters 


P 


DigZ3i] L,, = 0.520 L. = 0.237 C, = 1.140 C., = 0.485 C. = 1.140 


Number of Iterations — 58 Computer Run-time — 328 seconds 


Termination Caused by gi) - 41) Ee 


Graphical Comparison of Desired and Actual Frequency Response 


|Gain| a. 
(db) Desired 


Actual 





0.12 0.46 0.80 (Hz) 


Figure 9 


Data for Example 3 


a3 


4. Example 4 


Figure 10 contains the results of the optimization of the same 
circuit used in example 3, this time specifying a straight-line approxi- 
mation to the desired frequency response. Again, an approximate solution 
is the best that can be achieved. The final curve matches the desired 
reasonably well in general shape and bandwidth. 

5. Example 5 

Figure 11 contains the results of the optimization of a five- 
element band-pass filter where non-ideal inductors are considered. The 
desired frequency response is realizeable with the given circuit config- 
uration using ideal elements. When a resistance of R. = 0.1 is included 
in series with both inductors, the optimized response is slightly lower 
in gain than desired, but still is a good approximation as shown in the 
graphical comparison. The same circuit was optimized with ie = 0.01 and 
with R. = 0.5. The former gave virtually an exact match to the desired 
response; the latter produced a response that was generally 3.5 decibels 
below the flat portion of the desired response, and somewhat degraded in 
shape. 

6. Example 6 

Figure 12 contains the results of the optimization of an eleven- 
element low-pass filter whose desired frequency response is described by 
equation (4.1), with n = 10. The optimization produced a close fit to 
the desired response (maximum deviation from the desired response was 0.27 
decibels), but at the expense of relatively long computer run-time. This 
example required slightly over 32 minutes and 99 iterations to produce 
the final solution. After 60 iterations, the same general response had 


been found, with a maximum deviation of 0.42decibels (and an approximate 
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Circuit Configuration 





Starting Values of the Variable Parameters 


L, = 1.0 Ly = 1.0 L, = 1.0 C, = 1.2 Cy = 1.0 C, = 1.2 
Final Values of the Variable Parameters 

L, = 0.187 Lo = 0.492 L. = 0.316 C, = 1.160 C,, = 0.492 C. 

Number of Iterations — 110 Computer Run-time — 670 seconds 


Termination Caused by ) OR ep 


Graphical Comparison of Desired and Actual Frequency Response 


Desired — 
Actual 





0.12 0.46 0.89 — SHz) 


Figure 10 
Data for Example 4 
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Circuit Configuration 





Starting Values of the Variable Parameters 


L, = 1.0 Lueealiey C 


” = 1.0 C 


ae 


] 2 


Final Values of the Variable Parameters 


Lem 0...30) L, = 0.745 C, = 0.877 C, = 0.345 


1 2 1 2 


Number of iterations — 48 Computer Run-time — 18/7 seconds 
Termination caused by yoga eee 
Graphical Comparison of Desired and Actual Frequency Response 


|Gain| 
(db), 






Desired — 
Actual 


-10 


0.07 a Ll 


Figure 11 
Data for Example 5 
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Circuit Configuration 


1 2 3 5 ] 2 3 
C. =| 
Final Values of the Variable Parameters 
L, = 1.551 L, = 1.762 L, = 1.452 La = 1.043 L, = 
C, = 1.889 Co = 1.750 C, = 1.355 Cy = 0.948 C. = 


Number of Iterations — 99 


Termination caused by gi) ~ git) < 10 


12 


OF 
hes 


-6 


697 
374 


Computer Run-time — 1930 seconds 


Graphical Comparison of Desired and Actual Frequency Response 


Actual 





Figure 12 


Data for Example 6 
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Desired — 


run-time of nineteen minutes). This demonstrates the price paid to 
achieve a higher degree of "Closeness of fit" of the actual to the 


desired response. 


C. PROBLEM AREA 
In attempting to determine the maximum capability of the program, 

a complex sixteen-element low-pass filter was attempted. Results showed 
an unusually long computation time for one iteration (approximately 200 
seconds per iteration). Investigation disclosed that about eighty per- 
cent of the computational time was spent in the tree-finding subroutine 
of the modified Calahan section. As a result, after using 720 seconds 
of computer run-time and obtaining only 3.25 iterations, the problem was 
abandoned as requiring excessive computation time for a reasonable solu- 


tion. 
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V. A GUIDE TO THE USE OF THE OPTIMIZATION PROGRAM 


In this section the general characteristics of the optimization 
program are first discussed, followed by examples which are explained in 
detail and illustrate the use of the program. Finally, several guidelines 


are offered that result from experience with the program. 


A. GENERAL CHARACTERISTICS 
The optimization program consists of three main parts: 


1. The gradient-projection search section 


2. The modified Calahan program, which produces the frequency response 
of the given network. 


3. Subroutine GRAD which acts as the interface between sections one 
and two, and computes values for the gradient vector and the cost 
function. 

The program is dimensioned to accept fifteen variable network para- 
meters in a circuit with as many as 100 passive elements and 30 nodes. 
The circuit to be optimized is topologically described using nodes, 
branch element types, and their values. For the variable parameters, the 
specified value is the “starting point"; for elements that remain fixed 
the actual value is specified. Provisions are made for use of constraints 
and a cost-function weighting matrix, h. Termination criteria and the 
desired frequency response are supplied by the user. Program output 
includes: final parameter values; final values of the cost function and 
the norm of the gradient; the total number of iterations; printer plots 
of frequency response, phase, and delay for the final element values. 


Total storage space required for the computer program is 120K bytes. 


B. EXAMPLES 


Figure 13 is a general explanation of the required and the optional 
data cards necessary for use of the program. The specific format for each 
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I 1. NVAR — the number of elements to be varied 


2. NCON — if NCON is set to 0, Card Ila is 
used to designate common upper and lower 
12-3 4%5 bounds; if NCON is set to the number of 
constraint equations, Card group IIb is 
used, and Card IIa is omitted. 


3. KWIT — Maximum number of iterations desired 
4. EPS4 — ec of equation (2.3) 
5. WM(1) — hy in the weighting matrix, and 


7 re J 
6 7 
indicates that Card group VIII will be 
I a provided; if hy; is set equal to 0, equal 
2 weighting of ‘Ad frequency response points 
will be used by the program 
9 Ila 


6. The common lower constraint 
7. The common upper constraint 


IIb IIb 8. This block of cards designates all con- 
stants in the constraint matrix, if required 


III 9. The number of passive elements 
9 10 Vide 13 18 bb 10. The number of active elements 
II 11. The number of nodes 
12. The positive input node number 
ie) ae ee 13. The negative input node number 
14. The positive output node number 
IV 15. The negative output node number 


7 IV This block of cards describes the circuit 
20 topology. The first n cards must be the 
elements to be varied. 


V 16. The "from" node number 
17. The "to" node number 
2|__2e 18. The element type (R, L, C) 
VI 19. The element value 
V 20. The number of frequency intervals 
23 VI 21. The minimum frequency (Hz) 
22. The maximum frequency (Hz) 


VII VII 23. This block of cards describes the desired 
frequency response in decibels 


VIII 24. This block of cards is required only if 
VIII the user desires to specify the weighting 
matrix, and is necessary if item (5.)-is 
greater than 0 


Figure 13 
Data Cards Necessary for Use of the Optimization Program 
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card (or group of cards), is shown in the two examples that follow. It 
should be noted that two of the entries in Card I serve dual functions. 
1. If the second entry (NCON, which indicates the number of con- 
straint equations) is set equal to zero, upper and lower bounds for all 
variable elements are specified with Card Ila, satisfying the equation 
a< p,<b, 15 1,...,n where a is the common lower bound and b is the 
common upper bound; if a positive integer is entered, Card group IIb 
must be provided by the user, specifying the constants in the constraint 
equations, omitting Card IIa. (The formation of the constraint matrix 
is discussed in subparagraph B.2. of this chapter. ) 
2. If the fifth entry in Card I is set equal to zero, all frequency- 
response points will be weighted equally, and Card group VIII is omitted; 
if any positive number other than zero,is entered, the program stores this 


value as the first element of the weighting matrix, h Card group VIII 


1° 
(containing K-1 cards) representing Nos-++ shy must then be provided as 
the last section of the data deck. 

Two examples will be given using a four-element filter. The first 
example demonstrates use of the program in the simplest fashion, while the 
second illustrates the use of more complex constraint and weighting matri- 
ces. 

1. Example ] 

The circuit to be optimized is shown in Figure 14(a). The load 
resistor, R, is to be held constant; all other elements are to be opti- 
mized to match a given frequency response. A maximum of fifty iterations 
and ¢« = 0°” (see equation (2.3)) are the termination criteria. Al] 

' frequency points are to be equally weighted; eleven frequency points will 


be matched over the frequency interval (0.10,0.20). The three variable 


4) 





CARD NUMBER DATA FORMAT 


1 3 0 50 .000001 _ O. (314,2E15.8) 
2 0.1 4. (Zeon) 

3 4041 43 4 (7 (12,1X)) 

4 [same lil 

5 Zsa ml ..5 

6 ao ee (2(12,1X) AT, 1X, F 10.0) 
7 3 4R 1.5 

8 10 (13) 

9 0.10 0.20 (2E15.8) 

10 -0.005018 

(El Shs) 
20 -10.78127 
(b) 
Figure 14 


Data Cards for Example 1 
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parameters are to be constrained by the inequality O.1 <p, <4.0, 
i= 1,2,3. Figure 14(b) lists the required data cards. An explanation 


of the information on the cards follows. 


Card 1. In sequence: 

3 - The number of variables, n. 

0 - The number of constraints; in this case, 0 signifies 
that all variables are limited by like constraints, 
the values of which are specified on the next card. 

50 - The maximum number of iterations. The program wil] 
terminate with normal output after the 50th iteration, 
if no other termination criteria has been reached. 

0.000001 - This designates e. 

0. - Assignment of zero to this field indicates that all 
frequency-response points are to be weighted equally; 
that is, h. =|. 1 sale eo ook. 

Card 2. This designates the lower constraint as 0.1 and the upper 
constraint as 4.0 for all variables. 

Card 3. In sequence: 
4 - The total number of passive circuit elements 
0 - The total number of active circuit elements 

- The number of nodes 

- The positive input node number 

- The negative input node number 

- The positive output node number 


Sw fF — £ 


- The negative output node number 

Cards 4-7. These cards describe each branch element, as specified 
earlier. Note that the card for the non-variable load 
resistor follows the cards for the three variable elements. 

Card 8. The number of frequency intervals (the number of frequency 
points to be matched minus one) is entered. 

Card 9. The lower and upper frequency limits are entered. 

Cards 10-20. These cards list the desired frequency-response values 
at equal intervals, from the lowest frequency to the highest. 
There must be K points specified. The values must be in 
decibels. A maximum of 100 frequency points may be entered. 
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CARD NUMBER DATA FORMAT 


1 3 7 100 .000001 1. (314,2E15.8) 
2 1 

3 0 

4 0 

5 1 

: (E18.8) 
22 im 
23 7. 
29 | 0. 
30 4041 43 4 aye), 
31 Dae tt 2200 
32 2 ae eS 
33 0 3L 1.5 Jasssaans 
34 3 4R 1.0 
3 10 (13) 
36 ie 2 (2E15.8) 
37 - 005018 

he 
47 -10.78127 
48 2 
49 35 

(E15.8) 
57 " 


Figure 15 


Data Cards for Example 2 
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2. Example 2 
The same circuit is to be optimized, in this case using specific 
constraint and weighting matrices. The constraints are: 
0.1 < Py < ..5 
Urs < Po < 2.0 
LA < Pg = $.0 
2Po < Po 
These equations must be put into the form 
* 


Ap -B > 0 (5.1) 


~w 


Proper manipulation of the equations yields the seven constraint 


inequalities shown below in matrix form. 





The weighting matrix to be used is 


y=. 2. MR AT Oe ole OO ear. | 


This weighting puts greater emphasis on a close fit at the center portion 
of the frequency response curve. ‘The number of elements in the weighting 
matrix must correspond to the number of frequency points, K. Figure 15 
lists the required data cards. An explanation of the information on the 


cards foltows. 


* 
This notation means that each component of the vector Ap - B is 
non-negative. — “7 


Card |. In sequence: 


3 - The number of variables, n. 
7 - The number of constraints. This corresponds to 
the number of rows in A, and cannot exceed 32 
without re-dimensioning the program. 
100 - The maximum number of iterations. 
0.000001 - This designates e. 

1. - This assigns a positive value other than zero to 
the first element of the weighting matrix, and 
indicates that specified weighting matrix values 
will be provided. 


Cards 2-22. These cards list the values of each element of 
the constraint matrix A, sequentially by row. 


Cards 23-29. These cards list each value of the constraint 
matrix B. 


Cards 30-47. No change from example 1. 


Cards 48-57. These cards list the ten remaining values of 
the elements of h. 


Cy GUIDERINES 

Experience permits the suggestion of several guidelines that may be 
helpful in using the program. 

1. Termination Criteria 

Appropriate ruse for the maximum number of iterations and for e, 

depend primarily on the closeness of fit desired. Higher values for the 
maximum number of iterations (100 or higher), coupled with lower values 
OT ¢ (io7° or smaller) will ensure a closer fit, but at the expense of 
computer run-time. 


Additionally, the value for y in equation (2.4) 


, 3d(p) | 


Lad 








. 





ap 
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is designated within the main program as "ECON" and is set at a value of 
nm This termination criterion may be altered by the user if desired. 
2. Closeness of Fit at Termination 

Closeness of fit of the final frequency response to the desired 
response is somewhat dependent on how closely the circuit is capable of 
approximating the desired frequency response. If a close enough fit is 
not achieved, the following steps are recommended. 

Step 1. Impose more stringent termination criteria and re- 
initialize with the final values just obtained. (This step should be 
omitted if it is clear that the circuit configuration is not capable of 
coming close to the desired response. ) 

Step 2. Try different starting points for the variables if the 
termination value of the cost function is high. 

Step 3. Try a different network configuration. 

In addition, the value N in equation (2.2), the cost-function 
equation, has been set at N = 2 for least-squares curve fitting. This 
may be changed to N = 4 to penalize more heavily large deviations from 
the desired response. Statements 99 and 101 in subroutine GRAD both 
contain "**2" which should be changed to "**4" to accomplish this alter- 
ation. 

3. Number of Frequency Points 

Enough frequency points should be specified to completely describe 
the desired response. Doubling the number of frequency-response points 
does not double the computer run-time, but requires an increase of only 
about five to twenty percent in time (the larger percentage applies to 
small circuits, the smaller percentage to large circuits). This is 
because most of the run-time used is spent producing the network function 


(poles and zeros); once this is done, computation of the frequency 
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response is relatively fast. Failure to adequately describe the desired 
response may produce poor results. Experience with a limited number of 
examples indicates that an adequate description is provided when fre- 
quency points are provided to at least ten decibels below the desired 
flat response for low- and band-pass filters, and spaced over a close 
enough interval so that "corners" are well-defined. 
4. Computation Time 

Complex circuits (those with many trees) and any circuit with 
more than about eight elements require run-time in excess of ten minutes 
(on the IBM-360 digital computer). Investigation has shown that this 
time requirement is a function that varies directly with the number of 
trees in the circuit. Appropriate alteration of the Calahan subroutines 
to improve or alter the tree-finding method should provide a decrease in 
computation time necessary to produce a reasonable solution for large 


C1weuUlts. 
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VI. COMPARATIVE RESULTS OF GRADIENT-PROJECTION AND 
PATTERN-SEARCH METHODS 


Optimization of filter frequency response using pattern-search 
methods is the subject of a Naval Postgraduate School thesis being 
concurrently written by Captain James Lau, USMC. A comparison of the 
results of attempts to optimize the frequency response of three networks 
using each of the programs is shown in Figure 16. The first and second 
comparisons were optimized to a totally realizeable response, while the 
third was to a straight-line approximation. In all cases, the pattern- 
search optimization program proved faster and capable of producing more 


accurate results. 
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GRADIENT PROJECTION PATTERN SEARCH 
POTER the RESULTS RESULTS 


Eight-element band-pass 


Final cost-function value 0.564 x 1077 0.220 x 107° 

Computation time 328 seconds 184 seconds 
Five-element low-pass 

Final cost-function value 0.458 x 107° 0.204 x 107° 

Computation time 133 seconds 49 seconds 
Five-element low-pass 

Final cost-function value ese oor 

Computation time 132 seconds 50 seconds 


Figure 16 


Results of Three Network Optimizations, 
Comparing Gradient-Projection and Pattern-Search Methods 
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VII. CONCLUSIONS 


The purpose of this investigation was to determine the feasibility 
of using minimization techniques to cause a given network configuration 
to match a desired frequency response as closely as possible. This 
method was to use the computer's speed coupled to optimization techniques 
to provide the filter designer with a tool that quickly determines the 
network element values that most accurately cause the desired frequency 
response to be reproduced. Two existing programs capable of 

1. minimizing a given function, and 

2. producing a given network frequency response 
were adapted to produce a program capable of accomplishing the above- 
Stated purpose. 

The results presented in Chapter IV and the comparisons in Chapter 
VI lead to the following conclusions: 

1. It is feasible to obtain realistic approximations to a 
given frequency response, using a circuit capable of attaining that 
response. 

2. Satisfactory solutions using gradient-projection techniques 
can be obtained for filters with fewer than about twelve elements. An 
unreasonable amount of time may be required in attempting to optimize 
larger circuits. 

3. The pattern-search method for minimization is shown to be 
superior to the gradient projection method in both time and accuracy, for the 


examples compared. 
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Recommendations for logical extensions of this investigation fol low: 
1. Investigate a means to determine the gradient vector more 
accurately and more quickly. 
2. Alter the analysis program so that the tree-generating 


procedure is used only when necessary, rather than in each iteration. 
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APPENDIX A 
PROCEDURE FOR FINDING THE PARABOLIC MINIMUM USED IN THE 
FUNCTION-SEARCH INTERPOLATION METHOD 


1. The general equation of a parabola is: 
2 . 
Axo +Bx+C = y (A.1) 


The unknown coefficients A, B and C that will cause a parabola to pass 
through any three known points (not in a straight line) can be found 


by solving three simultaneous equations: 


2 


AX, + Bx, +C = yy (A.2a) 
2 _ 

Ax. + Bx, +€ = Yo (A.2b) 
2 a 

AX. + Bx. +C = Y3 (A.2c) 


2. These equations can be solved using the matrix equations: 


A X, X, ] yy 
B = Face cx 1 (Aes) 
ic 2 2 Y2 
Z 
c Xo Xo ] Y3 


3. The minimum of this parabola can be found by setting the first 
derivative of (A.1) with respect to x equal to zero and solving for ne 


d 2 a a 
aT ( Ax” + Bx +C-y) = 2Ax_ +2 See 
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APPENDIX B 
FUNCTION-SEARCH INTERPOLATION EXAMPLES 


1. Shown below are two examples (many others were noted, but are not 
included herein) where the gradient-projection minimization process 
seemed to be approaching a minimum, (no conclusive proof of this is 
offered, only the indications listed herein) when a distinct change in 
all of the variables as well as the cost function occurred. In the first 
four of the five iterations shown, the changes in the parameter values 
were small, and the cost function was being very slowly reduced. This 
behavior is characteristic of gradient methods when a minimum has been 
found (this can also be representative of moving down a "ridge line" or 
a "narrow valley"). This “jump" would not take place if the gradient 
method (as opposed to the gradient-projection method) were used. (Only 
five iterations are shown, as those before and after are not pertinent 
ToOmuhem ILLUS trabdena) 


First Example. Five-element (four-variable) low-pass filter 
Variable Parameter Values 


Iteration Cost 
number Py ip ne Pa function 
19 leant” 129541 -Orgo54" 0.5633 0.2425 
20 1.3284 1.9489 0.9624 0.5627 #£x0.2363 
2] Wass22e 1.9178 0.980), 0.5493 0.217] 

22 1.3468 1.9108 0.9743 0.5489 0.1978 
Ze wad SOmee] 65506 1.0565. 0.4367 02567 


Second Example. 


Variable Parameter Values 


Six-element (five-variable) low-pass filter 


Iteration Cost 

number Py Po BS P4 Ps function 
32 1.7525 1.8321 1.4150 1.0029 0.4440 0.02815 
33 1.7493 1.82708 1.4163 1,004) 024469 O2027c2 
34 1.7503 1.8269" 1.4185 oleO37 0.4488 0.02733 
25 1.7497 1.82638" 1.4182 “1.0039 0.450355 0-02720 
36 1.6674 1.6507 1.5231 1.0256 O7S@5see50 OMimec 
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2. These examples tend to substantiate the conclusion that the function- 
search interpolation method used in the gradient projection minimization 
program increases the likelihood of converging to the absolute minimum, 
for when the vector in the negative gradient direction intersects a 

point with a lower cost-function value — perhaps relatively distant from 
the base point — this new smaller-valued location will be taken up as 


the new base point. 
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